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ABSTRACT 

We have measured the redshift evolution of the density of Lyman hmit systems (LLS) in the in- 
tergalactic medium over the redshift range < z < 6. We have used two new quasar samples to 
(1) improve coverage at z ~ 1, with GALEX grism spectrograph observations of 50 quasars with 
0.8 < Zem < 1-3, and (2) extend coverage to z ~ 6, with Keck ESI spectra of 25 quasars with 
4.17 < Zem < 5.99. Using these samples together with published data, we find that the number density 
of LLS per unit redshift, n{z), can be well fit by a simple evolution of the form n{z) — n3.5[(l + z)/4.5]^ 
with ri3.5 — 2.80 ± 0.33 and 7 = 1.94^Qi32 ^r the entire range < z < 6. We have also reanalyzed 
the evolution of damped Lyman alpha systems (DLAs) in the redshift range 4 < z < 5 using our 
high-redshift quasar sample. We find a total of 17 DLAs and sub-DLAs, which we have analyzed in 
combination with published data. The DLAs with log HI column density > 20.3 show the same red- 
shift evolution as the LLS. When combined with previous results, our DLA sample is also consistent 
with a constant J^dla = 9 x lO""* from z = 2 to z = 5. We have used the LLS number density evolu- 
tion to compute the evolution in the mean free path of ionizing photons. We find a smooth evolution 
to z ~ 6, very similar in shape to that of Madau, Haardt & Rees (1999) but about a factor of two 
higher. Recent theoretical models roughly match to the z < 6 data but diverge from the measured 
power law at z > 6 in different ways, cautioning against extrapolating the fit to the mean free path 
outside the measured redshift range. 

Subject headings: cosmology: observations — early universe — quasars: absorption lines — galaxies: 
evolution — galaxies: formation 



1. INTRODUCTION 

Intergalactic absorption lines that are optically thick 
to radiation at the Lyman continuum edge are generally 
referred to as Lyman limit systems (LLS). The LLS are 
perhaps the single most valuable tracer of the evolution 
of the intergalactic medium since they provide a nearly 
direct measure of the mean free path of ionizing radiation 
at high redshifts. They are therefore a key ingredient 
in estimating the required sources of ionizing radiation 
needed to maintain the ionization in the intergalactic gas 
(Madau, Haardt & Rees, 1999 [MHR]). The expected 
rapid positive evolution in the number of such systems 
as the intergalactic gas becomes more opaque means that 
they could, in principle, be an indicator of the approach 
of the epoch of reionization. 

LLS are also the easiest systems to identify in quasar 
spectra since the abrupt break in the spectrum produced 
by the absorption systems can easily be seen in relatively 
noisy and low resolution observations. Since Tytler's pi- 
oneering work on the cosmological evolution of the LLS 
(Tytler 1982), numerous groups have measured LLS evo- 
lution over a range of redshifts extending from low red- 
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shift systems, measured with HST, to z ~ 5 (Bechtold 
et al. 2002; Sargent, Steidel & Boksenberg 1989 [SSB]; 
Lanzetta 1991; Storrie-Lombardi et al. 1994; Stengler- 
Larrea et al. 1995 [SL95]; Lanzetta et al. 1995; Peroux 
2001; Peroux et al. 2003 [P03]). The simplicity of iden- 
tifying LLS means that we can extend this type of mea- 
surement yet further, out to z ■^ 6, and this is the pri- 
mary goal of the present paper. Indeed LLS are the only 
neutral hydrogen line systems that we can readily trace 
beyond z '--^ 5. At these high redshifts, the lower column 
density Lya forest is too blended and the quasar contin- 
uum too ill-defined for these systems to be studied, and 
we can no longer determine the radiation damping wings 
that allow us to identify the stronger systems such as the 
damped Lya aborbers (DLAs). The LLS are therefore 
the one way we have to normalize the neutral hydrogen 
column density function at these redshifts, a crucial in- 
put to any modelling of the intergalactic medium. 

In the present paper we combine data on the LLS from 
a high redshift (z — 4 — 6) sample with a new low redshift 
sample (z < 1) and with pre-existing observations at in- 
termediate redshifts to compute the evolution of the LLS 
over the redshift interval z = — 6. The high redshift 
sample is based on high signal to noise obervations of a 
sample of 25 quasars (section 2.1). The new z < 1 sample 
is based on a large sample of z ^^ 1 quasars selected from 
grism spectra taken with the GALEX satellite. We use 
this to reinvestigate the low redshift density of LLS since 
the GALEX sample should be more homogeneous and 
unbiased than the previous HST and lUE samples used 
for this purpose. We find a lower value for the LLS den- 
sity at this redshift with the new sample. Remarkably we 
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TABLE 1 
High redshift quasars 



QSO 


^cm 


Exp (hr) 


^LLS 


'T'LLS 


•^^mctal 


2DLA 


log NuLA 


Rcf." 


FIRS0747+2739 


4.17000 


4.42 








3.4228 
3.9012 


20.80 
20.50 


20.85 (PW) 
20.50 (PW) 


PSS0209+0517 


4.17000 


3.00 


3.99060 


4.6 


3.98828 


3.6661 
3.8635 


20.48 
20.43 


20.45 (PW) 
20.55 (PW) 


PSS0926+3055 


4.19000 


3.00 














BRJ0426-2202 


4.32000 


3.00 














BR0334-1612 


4.36300 


3.00 


4.23037 


2.8 


4.23454 


3.556 


21.15 


21.0 (P03) 


PSS0747+4434 


4.39000 


3.00 


4.29520 


1.1 


4.27546 


3.762 
4.020 


20.15 
21.0 


20.3 (P03) 
20.6 (P03) 


BRI0952-0115 


4.42000 


4.25 


4.37745 


1.1 


4.37813 


4.0239 


20.60 


20.55 (PW) 


PC0953+4749 


4.46000 


2.00 


4.32691 


>4.2 


4.34002 


3.405 
3.889 
4.2444 


21.05 
21.03 
20.95 


21.15 (PW) 
21.20 (PW) 
20.90 (PW) 


BR1033-0327 


4.51000 


2.00 


4.18764 


>5.2 


4.17502 


4.174 


20.08 




BR2237-0607 


4.55000 


11.56 


4.25832 


>5.7 


4.24740 


4.080 


20.53 


20.52 (PW) 


PSS1347+4956 


4.56500 


1.20 


4.28753 


1.9 


4.27248 








BR0353-3820 


4.58000 


2.50 


4.36619 


1.5 


4.35953 








BR1202-0725 


4.61000 


3.00 


4.48567 


>6.2 


4.48012 


4.383 


20.62 


20.60 (PW) 


SDSS2200+0017 


4.78000 


6.83 


4.62166 


> 4.4 


4.65974 








SDSS0206+1216 


4.81000 


3.00 


4.34606 


>4.3 


4.32516 


4.3255 


20.78 




SDSS1737+5828 


4.85000 


9.00 


4.74902 


> 5.0 


4.74314 


4.7415 


20.70 


20.65 (PW) 


SDSS0211-0009 


4.90000 


6.00 


4.67670 


>4.1 


4.64357 


4.644 


20.18 




SDSS0338+0021 


5.01000 


7.25 


4.97626 


1.7 


4.97778 










SDSS1204-0021 


5.07000 


3.00 




















SDSS0231-0728 


5.42000 


4.33 


5.40306 


>4.6 


5.38874 














SDSS1044-0125 


5.74500 


5.75 


5.64344 


>2.6 
















SDSS0002+2550 


5.80000 


5.75 


5.72026 


>3.6 
















SDSS0836+0054 


5.82000 


11.83 


5.47923 


>2.5 
















SDSS141H-1217 


5.93000 


6.50 


5.89278 


2.3 


5.89602 














SDSS1306+0356 


5.99000 


9.75 


5.97264 


2.6 













Previously published values of logiVoLA. P03: Peroux et al. (2003); PW: Prochaska et al. (2007) 



find that the evolution over the full reshift range can be 
fit by a single power law of the form 2.80[(1 + z)/4.5]^'^^ 
and that this extends smoothly to near z = 6, beyond 
which the Lyman forest portions of the quasar spectra 
become too opaque to continue the measurement. 

We have also measured the DLAs (logA^ > 20.3) and 
sub-DLAs (here defined as 20.0 < logTV < 20.3) at z = 
4 — 5 in our high-redshift quasar sample, and used these 
in combination with previous work to test if the form 
of the neutral hydrogen column density distribution is 
changing as a function of redshift. The combined data 
are consistent with an invariant form of the function. 

We also compute the neutral hydrogen density in the 
range < z < 5 and find that J^dla is invariant over 
the redshift range 2 < z < 5. The LLS density evolution 
implies that this result may extend to z = 6 if the neutral 
hydrogen column density distribution function remains 
invariant. 

We use the form of the LLS evolution to compute the 
mean free path for ionizing photons, showing that the 
present observations decrease this compared with recent 
analyses (e.g. Faucher-Giguere et al. 2008). Finally we 
compare the evolution of the LLS and the mean free path 
to recent theoretical models. We assume Hq — 70, fim ~ 
0.3 and J^a = 0.7 throughout. 

2. DATA 

2.1. High redshift quasars 

We used a data set (Table 1) consisting of spectra 
of 25 quasars selected to span the redshift range from 
4.17 < Zem < 5.99. The sample was chosen from among 



the brightest quasars at these redshifts based only on ac- 
cessibility at the time of the observations. One additional 
BALQSO was observed in the program but has been 
eliminated from the present sample. While color-based 
quasar selection can introduce a bias toward quasars with 
strong LLS at lower redshifts, the La forest is sufficiently 
strong at these redshifts to ensure that the quasars would 
be picked out irrespective of the presence of an LLS (see 
Peroux 2001). The spectra were taken on the Keck II 
telescope with the ESI instrument in echellette mode 
with exposure times ranging from just under 2 hours 
to just under 12 hours. The resolution in this config- 
uration is comparatively low, ^ 5300 for the 0.75" slit 
width used, but the red sensitivity is high and the wave- 
length coverage is complete from 4000 A to 10, 000 A. 
At the red wavelengths, precision sky subtraction is re- 
quired since the sky lines can be more than two orders 
of magnitude brighter than the quasars. In order to 
deal with this issue, individual half-hour exposures were 
made, stepped along the slit, and the median was used 
to provide the primary sky subtraction. The frames were 
then registered, filtered to remove cosmic rays and arti- 
facts, and then added. At this point a second linear fit 
to the slit profile in the vicinity of the quasar was used 
to remove any small residual sky. This procedure gives 
an extremely accurate sky subtraction and a precise zero 
level in the spectra. The observations were made at the 
parallactic angle and flux calibrated using observations 
of white dwarf standards scattered through the night. 
These standards were also used to remove telluric ab- 
sorption features in the spectra, including the various 
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TABLE 2 

GALEX QUASARS 



QSO 



R.A. (J2000.0) 



Dccl. (J2000.0) 



^cni 



^LLS 



GALEX1713+6001 

GALEX1238+6143 

GALEX1438+3446 

GALEX1050+5819 

GALEX1421+5241 

GALEX1436+3539 

GALEX0040-4403 

GALEX0040-4340 

GALEX1434+3529 

GALEX1437+3508 

GALEX1237+6211 

GALEX1051+5838 

GALEX1434+3456 

GALEX0038-4352 

GALEX1052+5851 

GALEXlOOl+0200 

GALEX1417+5230 

GALEX0333-2817 

GALEX1723+5951 

GALEX1238+6202 

GALEX1423+5246 

GALEX1712+6005 

GALEX1237+6217 

GALEX0334-2743 

GALEX0332-2740 

GALEX1717+5947 

GALEX1420+5216 

GALEX1719+6013 

GALEX1436+3525 

GALEX1436+3453 

GALEX1724+5921 

GALEX1436+3537 

GALEX1240+6223 

GALEX1437+3442 

GALEX1435+3457 

GALEX1418+5223 

GALEX1234+6233 

GALEX0333-2725 

GALEX1435+3504 

GALEX1418+5252 

GALEX0334-2745 

GALEX1720+5932 

GALEX1002+0230 

GALEX1053+5855 

GALEX1054+5852 

GALEX1239+6206 

GALEX1712+6007 

GALEX1416+5259 

GALEX0040-4409 

GALEX1419+5232 



258.29000 
189.70300 
219.54100 
162.50100 
215.30099 
219.17599 

10.17800 

10.11710 
218.68699 
219.25599 
189.49800 
162.77200 
218.69799 

9.53954 
163.072998 
150.46299 
214.39999 

53.41400 
260.79299 
189.56700 
215.76699 
258.14498 
189.27900 

53.53450 

53.11040 
259.38699 
215.035003 
259.96398 
219.10299 
219.13200 
261.21398 
219.074005 
190.13900 
219.43200 
218.78199 
214.66000 
188.50500 

53.35220 
218.86799 
214.50700 

53.57030 
260.13501 
150.54499 
163.29200 
163.58700 
189.88099 
258.13198 
214.20100 

10.24390 
214.94299 



h 60.024200 
h 61.73180 
h 34.76800 
h 58.31790 
h 52.69649 
h 35.65800 
-44.061401 
-43.67160 
h 35.48809 
h 35.13760 
h 62.18399 
h 58.64239 
h 34.93420 
-43.87810 
h 58.85630 
h 2.0090301 
h 52.50840 
-28.29000 
h 59.85160 
h 62.0356 
h 52.77500 
h 60.091800 
h 62.28400 
-27.72710 
-27.67630 
h 59.79650 
h 52.27980 
h 60.22259 
h 35.42699 
h 34.89020 
h 59.35920 
h 35.62390 
h 62.39090 
h 34.71490 
h 34.95679 
h 52.39989 
h 62.55410 
-27.43070 
h 35.075801 
h 52.86690 
-27.75120 
h 59.53530 
h 2.50739 
h 58.92580 
h 58.87179 
h 62.10520 
h 60.12210 
h 52.98410 
-44.16410 
h 52.54639 



0.8179 
0.8323 
0.8331 
0.8350 
0.8438 
0.8450 
0.8524 
0.8726 
0.8729 
0.8800 
0.9129 
0.9150 
0.9269 
0.9359 
0.9700 
0.9705 
0.9877 
0.9906 
0.9964 
0.9976 
0.9993 
0.9993 
1.0190 
1.0281 
1.0396 
1.0482 
1.0600 
1.0600 
1.0648 
1.0676 
1.0770 
1.0776 
1.1115 
1.1146 
1.1174 
1.1202 
1.1231 
1.1400 
1.1472 
1.1600 
1.1634 
1.1691 
1.1691 
1.1800 
1.1835 
1.2008 
1.2037 
1.2325 
1.2325 
1.2469 



0.7784 



0.8864 



0.7935 



0.7770 
0.8428 



1.1100 

0.8325 
0.7651 

0.8201 



atmospheric bands. The final extractions of the spectra 
were made using a weighting provided by the profile of 
the white dwarf standards. Wavelength calibrations were 
obtained using third-order fits to CuAr and Xe lamp ob- 
servations at the begining and end of each night, and 
the absolute wavelength scale was then obtained by reg- 
istering these wavelength solutions to the night sky lines. 
The wavelengths and redshifts are given in the vacuum 
heliocentric frame. A sample of the final spectra of the 
25 quasars is shown in Figure 1 and all of the spectra are 
given in the Appendix. 

Each spectrum was visually inspected to determine the 
wavelength (Alls) below which the spectrum dropped by 
more than a factor of 2.7 (r = 1). We adopt the redshift 
-^LLS = Alls/912 — 1 as the Lyman limit redshift. The 
procedure used is shown in Figure 2 (see also Appendix) 
where we plot the quasar spectra in the vicinity of the 



LLS and show the average fluxes in wavelength intervals 
above and below the position of the LLS. The results 
are summarized in Table 1, where we give the quasar 
name, its emission redshift Zem, the ESI exposure time 
and the LLS redshift (columns 1-4). Where no LLS 
is detected to our adopted wavelength limit of 4500A 
{ziim — 3.923) we show a blank entry in Table 1. For 
each Lyman limit system we searched for corresponding 
metal-line absorption and where this was detected we 
also give the metal- line redshift (column 6). The optical 
depth was now measured by computing the actual flux 
for a rest frame interval of 30 A immediately below the 
LLS wavelength and taking the log of the ratio of this 
quantity to the flux obtained by extrapolating a power- 
law fit to the quasar spectrum at 10 ~ 160 A longward 
of the LLS, excluding the region around the O VI -I- 
Ly/3 quasar emission. The procedure is illustrated in 
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Fig. 1. — ESI spectra of the QSOs used in this investigation. 



THE EVOLUTION OF LYMAN LIMIT ABSORPTION SYSTEMS TO REDSHIFT SIX 



2.0x10"'' 
1.5x10"" 
1.0x10"" 
5.0x10"" 

-5.0x10"" 





L PSS0209+0517 i 


: REDSHIFT=4.17 | || 


II , i 1 


: z^=3.990 1, 

^" 1' 11 ill , 


1,1 ' 


iJ 


- 


i,;iii,i!ili il 




wy, 


vrfV 




1 


iij 


ii^ 


k 


'|i' III III 







800 850 900 950 1000 

VACUUM HELIOCENTRIC WAVELENGTH (ANGSTROM) 




800 850 900 950 1000 

VACUUM HELIOCENTRIC WAVELENGTH (ANGSTROM) 



Fig. 2. — Lyman limit systems with r > 1 in quasars at inter- 
mediate redshift. The spectra are plotted in the rest frame of the 
quasar. The blue vertical line marks the wavelength Alls of the 
LLS and the redshift zlls is marked alongside. The red horizontal 
bars show the average flux in 30 A bins below Alls S'ld to redward 
of -^LLS + 10 A. The blue dashed line shows the predicted contin- 
uum level below Alls based on a power-law fit to the redward 
continuum. 

Figure 2 and the optical depths are tabulated in column 
5 of Table 1. 

We tested the accuracy of the LLS redshift determina- 
tions and optical depth measurements by simulating LLS 
in the spectra and then measuring the recovered redshift 
and optical depth. For each quasar for which there was 
sufficient wavelength space between the first real LLS 
and the quasar La emission line, we inserted a false LLS 
at a randomly generated wavelength and with a specified 
optical depth. Below the wavelength of the false LLS we 
added in artificial noise matched to the quasar exposure 
time. We then remeasured the LLS redshift and optical 
depth using our standard procedure. We repeated this 
twenty times for each quasar for three optical depths: 
r = 1.0, 1.5 and 2.0. 

The measured redshift offset as a function of LLS red- 
shift is shown in Figure 3 for t = 1.0 and 2.0. In both 
cases we have shown the individual offsets of each simula- 
tion with smaller symbols, and the average in a particular 
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Fig. 3. — Redshift offset of simulated Lyman limit systems as a 
function of redshift for simulated LLS optical depth, r = 1.0 (top) 
and 2.0 {bottom). Black squares show individual test points for 
the quasars in our high redshift sample. Large red squares are the 
mean for each quasar. 



quasar with the larger symbol. As would be expected, 
offsets tend to be positive, i.e. the measured redshift is 
very slightly higher than the real one. This is a conse- 
quence of blending with the forest. The effect is larger 
for lower optical depth and at higher redshift, where the 
forest becomes richer. Typical offsets in the r = 1.0 case 
where the effect is largest are ~ 0.01 below z — 5, rising 
to ~ 0.02 - 0.03 in the highest redshift quasars. The er- 
rors are generally smaller for the r > 1 systems. We shall 
consider the systematic errors which can be introduced 
by these offsets in the next section. However, the sys- 
tematic errors introduced by the redshift uncertainty are 
smaller than the statistical uncertainty. We also used the 
simulation to calibrate the measurements of the optical 
depth and examine the systematic errors in the deter- 
mination of this quantity. We find that the errors are 
dominated by the continuum fitting, with a dispersion of 
0.3 in T, and that there is a systematic offset of +0.3 in 
T. This is introduced by the use of a single power law fit 
which systematically overestimates the continuum. We 
have corrected the values of r given in Table 1 to allow 
for this offset. The statistical errors and other system- 
atic errors are completely dominated by the continuum 
fitting error. We consider the effects of the continuum 
fitting error on the determination of the LLS number 
density in the next section. 

Finally we searched for all damped and sub-damped 
Lyman alpha systems in each spectrum and give the 
Lyman-alpha redshift and column density for each de- 
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tected system with iVn > 10 cni~ (columns 7 and 8). 
Where more than one DLA was detected we give the 
properties of ah the lines. Where the DLA has a pre- 
viously measured value in the literature we give this in 
column (9) together with the reference. 



2.2. 



1 GALEX quasars 



To improve the determination of the LLS density at 
low redshift we selected a sample ofz = 0.8 — 1.3 quasars 
from the table of Cowie et al. (2010) who tabulated all 
of the Lyman alpha emission-line objects in the 9 deep- 
est GALEX blank-field grism spectrograph observations 
covering 8.2 square degrees of sky. The GALEX spectra 
were obtained from the Multimission Archive at STScI 
(MAST). Morrissey et al. (2007) detail the spectral ex- 
traction techniques used by the GALEX team in extract- 
ing the spectra from the grism data. We restricted our- 
selves to quasars lying at redshifts z = 0.8 — 1.3 and with 
fluxes above 8 x 10~^^ erg cm~^ s~^ A~^ where the Ly- 
man breaks can easily be picked out. The sample is NUV 
(1900-2800 A) magnitude selected (see Cowie et al. 2010 
for a detailed description of the selection function) and 
we do not expect any significant biasing from the pres- 
ence of a LLS at z < 1.3 (2090 A) which would only 
minimally reduce the NUV magnitude. The GALEX 
sample used is given in Table 2, where we give the name 
of the object, the quasar redshift Zem, its RA(2000) and 
Dec(2000) and the redshift of the LLS if one is detected. 
The observed wavelength range means that only systems 
with zlls > 0.54 are included. 

We inspected all of the GALEX spectra for these ob- 
jects and picked out those with strong LLS (r > 1). To 
do this we fitted a power law continuum to the spec- 
trum above 2100 A and then found the first location at 
which the spectrum dropped by more than a factor of 2.7 
and did not subsequently recover at shorter wavelengths. 
The GALEX spectra in which we detected such a strong 
LLS are shown in Figure 4 (see also Appendix), where 
we order the sample by the red continuum flux. An im- 
portant point to note is that the GALEX spectra consist 
of two wavelength segments: one corresponding to the 
FUV grism (~ 1400 - 1800 A) and one to the NUV 
grism (- 1900 - 2800 A). We can stfll select LLS in the 
z = 0.97 — 1.08 range since the flux break across the gap 
is easily seen; consequently, the search pathlength is the 
redshift range including the gap. In the case where a 
LLS did fall in the gap, we would not be able to make a 
precise measurement of the break redshift and we would 
have to address the redshift uncertainty in the analy- 
sis. In fact we do not find any LLS in our sample whose 
redshift lies in the gap. LLS systems were found in 9 
of the 50 quasars at redshifts above the z = 0.54 limit 
where we could have detected such a system, though one 
system (in GALEX 1418-1-5223) lies close to this region 
and we consider the changes in the results that would 
be produced if it lies at a lower redshift. The LLS in 
GALEX 1724-1-5921 is somewhat uncertain because of 
the noisy spectrum. We examine the reduction in the 
number density which would be produced by removing 
it in the analysis. 

3. LYMAN LIMIT SYSTEMS 
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Fig. 4. — GALEX quasars with Lyman limit systems. The spec- 
trum of the quasar is shown in the rest frame with the position 
of the quasar emission lines marked in red. The short wavelength 
portion of the spectrum is from the GALEX FUV grism and the 
long wavelength portion from the GALEX NUV grism. The posi- 
tion of the LLS is marked in purple. The figures are ordered by 
the red continuum flux above the Lya line. 



TABLE 3 
Parameters of Figure 5 



Redshift bin 


#QSO 


#LLS 


n{z) 


Error 


(z) 


(0.36,1.1)" 


51 


11 


0.63 


0.19 


0.64 


(0.58,1.25)'' 


49 


7 


0.30 


0.12 


0.82 


(0.36,1.25)'= 


100 


18 


0.44 


0.10 


0.72 


(2.4,2.9) 


77 


29 


1.84 


0.33 


2.74 


(2.9,3.8) 


106 


52 


2.43 


0.33 


3.20 


(3.8,4.2) 


52 


31 


3.77 


0.64 


4.04 


(4.2,5.0) 


61 


33 


4.04 


0.70 


4.39 


(5.0,6.0) 


8 


6 


8.91 


3.49 


5.69 



^ Stcngler-Larrea et al. 1995 sample 

^ GALEX sample 

^ Combined Stengler-Larrea and GALEX samples 

We first computed the number density of LLS per unit 
redshift path dz, which we write as n{z). As is con- 
ventional, we excluded LLS lying within 3000 km s~^ of 
the quasar and computed the redshift path accordingly. 
The data can be analysed assuming a constant density 
of sources in a redshift interval or as a fit to an assumed 
redshift evolution over the entire redshift interval. To 
determine n{z) and its 1 a error in each redshift interval 
we used the maximum likelihood fit (MLF) in the form 
given by SSB. The corresponding results are shown in 
Figure 5. For the higher redshift interval we have shown 
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only the full samples (the present data and all previous 
data taken from the compilation of P03). However, for 
z < 1 we have shown the results obtained separately from 
the present GALEX sample (blue diamond) and the older 
SL95 sample (red triangle) as well as that for the com- 
bined sample. We have chosen to show the independent 
values since it is possible that the difference could reflect 
systematic selection effects in the archival UV observed 
quasar samples analysed by SL95 which might result in 
an upward bias in the LLS density. However, the SL95 
and GALEX samples are consistent within the statistical 
uncertainties and we shall use the combined value in the 
following text. The present data roughly doubles the low 
redshift sample and provides a corresponding reduction 
in the error over the previous SL95 result (Table 3). The 
values of n{z) in each redshift interval are summarized 
in Table 3, where we give the redshift bin, the number 
of QSOs used, the number of LLS detected, the value 
of n{z) and its 1 a error, and the average redshift of 
the LLS. For the combined sample, we find a value of 
n{z) = 0.46 ± 0.11. If the LLS in GALEX 1724+5921 
is eliminated the values in the z = 0.36 — 1.25 interval 
would reduce to 0.44 ± 0.10, while if the redshift of the 
LLS in GALEX 1418+5223 is lowered to 0.97, the value 
is unchanged at the quoted level of significance. We have 
eliminated the LLS in GALEX 1724+5921 in the subse- 
quent analysis and the numbers in Table 3 correspond to 
the sample excluding this quasar. 




2 3 4 5 6 7 

1 +z(LLS) 



Fig. 5. — Filled squares: Density of Lyman limit systems as a 
function of mean Lyman limit system redshift for the complete 
sample of quasars. Error bars are 1 a. For each bin, the number 
of quasars, number of LLS, the density of LLS and the mean LLS 
redshift are given in Table 3. Maximum likelihood fits are shown for 
1.5 < z < 6 (green solid line) and < z < 6 (green dashed line). At 
low redshifts we show the values inferred from the GALEX sample 
alone (blue diamond) and the SL95 sample alone (red triangle) as 
well as the value from the combined data sets, which is shown by 
the black square. The black dashed line is the fit of SL95. 

We also examined the effects of the estimated redshift 
and optical depth errors derived from the simulation and 
discussed in the previous section and shown in Figure 3. 
For z = 4.2 - 5, reducing the LLS redshift by the offset 
—0.01 for the t = 1 case would reduce the measured n(z) 
from 4.04 ± 0.70 to 3.77 and at z = 5 - 6, reducing the 
LLS redshift by -0.03 would change n(z) from 8.91+3.49 
to 7.03. The effect is smaller for higher optical depth 



systems and the systematic changes are smaller than the 
statistical uncertainty. The 0.3 scatter in optical depth 
introduces an Eddington bias since there are more low-r 
systems to scatter up than scatter down. However, this 
effect is also very small, a few percent in n{z), compared 
with the statistical error. 

In analysing the redshift evolution of the LLS it is usual 
to parameterize the evolution in n(z) as 



n{z) = no(l + zy 



(1) 



where ng is the value of n(z) at z = 0, and to estimate 
the parameters using a maximum-likelihood fit to the 
full data set in a given redshift interval. The motivation 
for this form lies in the historical description, based on 
assuming that LLS lie in galactic haloes, but it is also 
a simple parameterization which provides a reasonable 
choice of form fitting for the sample and is useful for 
comparing with previously derived results. However, the 
use of uq results in a strong coupling between this pa- 
rameter and 7, reflecting the lever arm in the power law. 
This may be avoided, and the display of the parameter 
fit simplified, by normalizing nearer to the center of the 
redshift range and we use the form 



n{z) = 713.5 



1 



4.5 



(2) 



where 713.5 is the value of n{z) at z = 3.5, to show the 
error contours on the parameters. 
We calculated the log-likelihood function 

lnL = a:lnno + 7^Ml + 4Ls)- (t^) 

xE[(l + ^^„X+^-(l + 4,J^+^] (3) 
i=i 

where the notation follows that of Peroux (2001) and x 
is the number of LLS and y is the number of quasars, 
and used this to estimate 7 and 713.5 a-nd their errors in 
the redshift intervals 1.5 < z < 6 and < z < 6. The 
results are summarized in Table 4, where we also list the 
parameters and fitting ranges of previous fits from the 
literature. 

The fits are shown in Figure 5 as the solid (1.5 < z < 6) 
and dashed (0 < z < 6) green lines, compared with the 
redshift evolution of the density of LLS in the complete 
quasar sample (filled black squares) in suitable redshift 
bins. In contrast with previous results, and primarily 
as a consequence of the new low-redshift 7i(z) value, the 
fits over the two redshift ranges are practically indistin- 
guishable and give consistent values for 7 and 713.5, as 
is seen in Figure 6, which plots the 1, 2, and 3 a error 
contours. A single parameterizaton (71,3.5 = 2.80 + 0.33, 
7 — 1.94to;32) fits well over the entire redshift range, 
< z < 6. The evolution inferred with this data set is 
generally steeper than previously measured (see Table 4), 
except for the result of Peroux (2001). In particular, 
the evolution is steeper than that measured by SL95 for 
0.32 < z < 4.11, (the black dashed line in Figure 5), as 
a result of the improved low-redshift determination from 
the large combined sample. 
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TABLE 4 
Parameters of Maximum Likelihood Fit 



Sample'^ 



Rodshift range 



QSO 



LLS*" 



7 



"3.5 



no 



This 


paper 


Rcf. 


(1) 


Ref. 


(2) 


Ref. 


(3) 


Ref. 


(4) 


Ref. 


(5) 



1.5 -6 
0.0-6 
0.40 - 4.93 
0.32 - 4.11 
0.40 - 4.69 
0.67 - 3.58 
3.3- 4.4 



216 
315 
148 
169 
90 
90 
469 



157 
175 
84 
80 
48 
54 
192 



2 04+0-62 

1 Q4+0-36 
^•^^-0.32 

1.50 ±0.39 

1.55 ±0.45 

0.68 ±0.54 

5.2 ±1.5 



2.84 
2.80 



+0.29 
0.37 

+0.33 
0.33 



0.13 
0.15 

0.06 
0.25 
0.27 



+0.15 
0.05 

+0.17 
0.10 

+0.20 
0.13 



0.76 



1.5 ±0.2 



References: (1) Poroux (2001); (2) Stcnglcr-Larrca ct al. (1995); (3) Storrio-Lombardi et al. (1994); (4) SSB (1989); (5) Prochaska ct al (2010). 



Refers to number of LLS with r > 1, except for Ref. 4 (r > 1.5) and Ref. 5 (r > 2). 
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Fig. 6. — 1 it, 2 c, and 3 a maximum likelihood contours of 
"3.5i 7 for LLS in the redshift intervals < zi^i^g < 6 {black filled 
square, contours) and 1.5 < zi^i^g < 6 {red triangle, contours). 



4. DAMPED LYMAN ALPHA SYSTEMS 

We searched the spectra for absorption Hnes with 
damping wings and found a total of 17 DLAs 
(log7V(HI) > 20.3) and sub-DLAs (20.0 < log7V(HI) < 
20.3) within our z < 5.1 quasar sample, 14 DLAs and 
3 subDLAs. We believe the sample is complete to the 
level logAf(HI) = 20. However, in the highest redshift 
quasars the lower column density systems are hard to 
measure and we have not attempted to extend the mea- 
surements below logA''(HI) = 20. DLAs in general be- 
come too hard to measure beyond z — b.l because of the 
thickness of the forest. We measured Af(HI) for these 
systems by fitting to the damping wings (Figure 7 and 
Appendix), and column 8 of Table 1 compares our mea- 
sured values of logA^(HI) with previous measurements 
from Prochaska et al. (2007) and P03. Fourteen of the 
systems have existing measurements. All of the DLAs in 
the present sample are shown in Figure 7 together with 
the damped Lq fits (red solid lines). Nine of the DLAs 
overlap with DLAs in the P03 sample and the P03 fits 
are shown with blue dotted lines in Figure 7. One DLA 
in the P03 sample is classified as a sub-DLA here (the 
z = 3.762 system in PSS0741+4434), and one DLA in the 
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Fig. 7. — Damped Lyman alpha systems in the high redshift 
sample. The red dotted line shows the fit corresponding to log Nn 
from Table 1. Where present, the blue dotted line shows the fit 
from P03 (see Table 5). 



P03 sample (the z = 4.06 system in SDSS0338-f002) is 
not confirmed. However, in general the measured column 
densities agree well with previous measurements. The 
present sample substantially increases the path length 
above z = 4.5. We summarize the present results com- 
bined with z > 4.5 quasars from P03 and Guimaraes et 
al. (2009) in Table 5. We have used these combined sam- 
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pies to compute the redshift evolution of DLAs out to 
z = 5.1 to compare with that of the LLS. 



1.00 



z 



0.10 



0.01 




DLAs, at least for z > 1. The situation is less clear for 
the higher column density (logA^(III) > 21) DLAs. At 
low redshifts {z < 3.5) the LLS can be scaled down by 
a factor of 33 to match the high column density DLAs 
but the errors at high redshift could permit a substantial 
deviation from this form. We adopt a single power-law 
fit, (fN{m)/dNdX ex N{m)-^, to the column den- 
sity distribution as a function of redshift between the 
LLS (log TV (HI) = 17.2) and DLA (log TV (HI) = 20.3 and 
21.0) column densities. For log A^(HI) = 20.3 the individ- 
ual points are consistent with a single index of /3 = 1.28 
for 1 < z < 5, whereas an index of /3 = 1.40 is required 
for logA''(HI) = 21.0. The steeper index over the wider 
column density range is a consequence of the well known 
turnover in the column density distribution function at 
high column density (e.g. Petitjean et al. 1993; Storrie- 
Lombardi & Wolfe 2000; Prochaska et al. 2005; O'Meara 
et al. 2007; Guimaraes et al. 2009). 



Fig. 8. — Number density of DLAs as a function of redshift for 
various samples and computed for two values of the limiting col- 
umn density, A'^h- For Ajj > 20.3, colored symbols denote: Rao 
et al. (2006; black open squares); Prochaska et al. (2005; red filled 
diamonds); Peroux ct al. (2005; green filled triangles); Guimaraes 
et al. (2009; blue filled triangles). For Nn > 21.0, the green open 
circles show values computed from Peroux et al. 2005. The filled 
black symbols show the high redshift sample of this paper (Table 
1) combined with all previous high redshift samples for Ajj > 20.3 
(squares) and Ajj > 21.0 (circles). The solid black line is the max- 
imum likelihood fit to the LLS with < z < 6 from the combined 
sample (see Table 4). The dashed and dot-dashed black lines are 
this relationship scaled by a factor of 7.5 and 33, respectively. 

In Figure 8 we compare the redshift evolution of the 
DLAs with those of the LLS. The black line shows the 
maximum likelihood fit to the evolution of chi/dX in the 
LLS over the range < z < 6, where dX is defined as 



-1/2 



^ = (1 + z)2 [(1 + ^)2(i + ^Oj^) _ z{2 + z)nA] 

dz '- ■' 

(4) 
and dn is the number of t > 1 LLS in the interval dX. 
The black squares show the observed number of DLAs 
with logiV(HI) > 20.3 in the redshift intervals [3.5,4.3], 
and [4.3,5.1], computed from the present data set com- 
bined with preceding measurements. The error bars are 
±1 a based on the Poisson errors corresponding to the 
number of systems in the bin (Gehrels 1986). We have 
also computed the errors using a jackknife estimation. 
These agree well for bins with higher numbers of sys- 
tems but underestimate the error in bins with only a 
small number of DLAs. We also show lower redshift 
DLA number density measurements from the literature 
(Prochaska et al. 2005, Peroux et al. 2005, Rao et al. 
2006 and Guimaraes et al. 2009.) The black circles show 
the distribution for the log A^ (HI) > 21 systems with the 
P03 results shown with green open circles. 

If the neutral hydrogen column density distribution 
function is invariant with redshift the redshift evolution 
of the LLS and the DLAs would have the same form in 
Figure 8. This appears to be the case for the log A^(HI) > 
20.3 systems where a simple scaling down of the LLS 
numbers by a factor of 7.5 replicates the evolution of the 
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Fig. 9. — Hdla as a function of redshift. Symbols have the same 
meaning as in Figure 8. 



Finally, we have computed JIdla for the combined 
sample for 3.7 < z < 4.4 and 4.4 < z < 5.1 and com- 
pare these values with results from the literature in Fig- 
ure 9 {n here includes both hydrogen and helium). The 
present data (black squares and 1 a error bars) show 
no evolution in J7dla to z = 5 (see also Prochaska et 
al. 2005; Peroux et al. 2005; Noterdaeme et al. 2009; 
Guimaraes et al. 2009). This remains true if we expand 
the sample to include the subDLAs, which raises J^dla 
by about 30% in the highest redshift bin. The evolution 
is consistent with a constant fl over the observed redshift 
range. 

5. DISCUSSION 

The density evolution of the LLS is most closely re- 
lated to the attenuation length of ionizing photons in 
the intergalactic medium since it is systems with column 
densities near the value of the LLS that dominate the 
neutral hydrogen opacity. Thus the conversion of e, the 
source function of ionizing photons, to F, the ionization 
rate in the intergalactic medium, can be determined if 
we know the evolution of the LLS and the shape of the 
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TABLE 5 

DLAS IN 2; > 4.5 QUASARS 



Quasar 


^em 


^min 


ZDLA 


logAr(H) 


Source" 


Other'' 


PSS1140+6205 


4.500 


3.113 


0.000 


0.00 


1 




PSS0808+5215 


4.510 


3.181 


3.114 


20.60 


1 




BR1033-0327 


4.510 


2.888 


4.174 


20.08 







PSS1723+2243 


4.515 


3.062 


3.697 


20.35 


1 




PSS1715+3809 


4.520 


3.004 


3.341 


21.05 


1 




BR0019--1522 


4.528 


2.970 


3.437 


20.92 


2 




PSS0134+3307 


4.536 


2.976 


3.760 


20.60 


1 


P 20.6 


BR2237-0607 


4.550 


2.941 


4.080 


20.53 





P 20.50 


PSS1347+4956 


4.565 


2.963 


0.000 


0.00 







BR0353-3820 


4.580 


3.022 


0.000 


0.00 







BR1202-0725 


4.610 


3.112 


4.383 


20.62 





P 20.49 


SDSS0210-001 


4.700 


3.177 


0.000 


0.00 


1 


P 20.5 (3.420)* 


BRJ0307-4945 


4.728 


3.138 


4.460 


20.80 


2 




SDSS2200+001 


4.780 


3.214 


0.000 


0.00 







SDSS 0206+121 


4.810 


3.007 


4.326 


20.78 







SDSS0941+594 


4.820 


3.322 


0.000 


0.00 


1 




PSS1247+3406 


4.897 


3.314 


0.000 


0.00 


1 




SDSS0211-000 


4.900 


3.255 


4.644 


20.18 







SDSS1737+582 


4.850 


3.309 


4.741 


20.70 





PW 20.65 


SDSS0338+002 


5.010 


3.480 


0.000 


0.00 





P 20.4 (4.060)* 


SDSS1204-002 


5.070 


3.437 


0.000 


0.00 







SDSS0756+410 


5.090 


3.551 


4.360 


20.15 


1 





'^ Source for cols. 1-5: 0: Present work; 1: Guimaraes et al. (2009); 2: Peroux et al. (2003) 
^ Other published values of log Ar(H) (redshift): P: Peroux et al. (2003); PW: Prochaska et al. (2007) 
Not confirmed. 



column density distribiution in the neighborhood of the 
LLS column density. 

Following the fundamental work of Madau, Haardt and 
Rees (1999; MHR), there have been a number of revis- 
its to the problem, updating the calculation using recent 
values of the evolution parameters and the column den- 
sity shape. The most recent version is given by Faucher- 
Giguere et al. (2008). At z > 2 these calculations are 
enormously simplified by the short mean free path of the 
photons which means the calculation is essentially local 
(MHR). In this limit, if we can approximate the column 
density distribution function in the neighborhood of the 
LLS column density as N~^ , the mean free path is simply 
related to the LLS system density n{z) as 



AI^vq, z)n{z) oc 



J a 



/o r-^(l~e-) 



or 



Al{iyo,z) 



,1-/3 



r(2-/3) H{z){l + z)n{z) 



(5) 



(6) 



where r is the optical depth at the Lyman edge, a is the 
optical depth above which n{z) is measured (here a ~ 1) 
and vq is the frequency of the Lyman continuum edge 
(see for example the equations in Miralda-Escude (2003) 
or Choudhury (2009) , and note that F in this equation is 
the gamma function and not the ionization rate). MHR 
adopted redshift-invariant values of /? = 1.5 and 7 = 2 
to calculate the values, basing these parameters on the 
work of Petitjean et al. (1993), Hu et al. (1995), and 
Kim et al. (1997) for the /? index. Faucher-Giguere et al. 
(2008) updated the parameter /3 to 1.39 using the results 
of Misawa et al. (2007), used Stengler-Larrea's fit to the 
LLS population, and again assumed a redshift-invariant 
set of parameters. 



The present work supports the assumption of redshift 
invariance, and our power law index fit of 0.28 to the cu- 
mulative distribution between the LLS and the log Nn > 
20.3 DLA column densities is not substantially different 
(in the differential form of 1.28) from the Misawa et al. 
(2007) value adopted by Faucher-Giguere. However, the 
present data implies a higher value of the LLS density 
at z > 2 than was adopted by Faucher-Giguere, with a 
corresponding decrease in the mean free path since, as 
was already known from the work of P03, the SL95 fit to 
the high redshift LLS is low because their sample did not 
extend to higher redshifts. Combining the evolutionary 
form of the LLS with equation (6), and using the high 

redshift approximation, H{z) = Hoflm (1 + z)^/^, gives 



Al{iyo,z) 



7I-/3 



(4.5) 



-2.5 



r(2 - 13) Hon 



1/2 



^3.5 



4.5 



-(2.5+7) 



(7) 



where we have used the normalization at z — 3.5 which 
places the errors primarily in the 7 exponent (see Fig- 
ure 6). For /3 — 1.3, the maximum likelihood fit for 
< z < 6 gives 



Al{iyo,z) = 50 



l + z 



4.5 



-4.44^ 



Mpc. 



(8) 



For /3 — 1.5, the normalization becomes 37 Mpc and for 
13^ 1.1, 61 Mpc. 

We show the mean free path of equation (7), and also 
that computed from the measured tt-lls at that red- 
shift, in Figure 10, where we compare with the mean 
values computed by MHR and Faucher-Giguere et al. 
(2008) and with the mfp determination of Prochaska 
et al. (2009a) in the range 3.6 < z < 4.2. Note 
that the 1 a errors shown are formal statistical errors 
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Fig. 10. — (Top) Mean free path as a function of redshift. The 
filled squares show the data points computed from Equation (7) 
using the maximum likelihood estimate of tills in the redshift 
bins indicated by the horizontal lines. 1 a errors are shown in the 
y direction, the points are computed for /3 = 1.28 but may be 
simply scaled to other choices of /3 using the functional form of 
Equation (7). The thick red line is the fit to our combined sample, 
from Equation (7) with 7 = 1.94, /3 = 1.28; the thin red lines show 
the 1 (T range in /3 (1.1 to 1.5). The black dot-dash line is the result 
of Faucher-Giguere et al. (2008), with 7 = 1.5, /3 = 1.39. We note 
that, in comparing with MHR (blue dotted line), Faucher-Giguere 
et al. did not correct for the cosmology, making the difference 
between their results and MHR appear larger in their Figure 4 than 
in the present figure. The green dashed lines are the power law fit 
of Prochaska et al. (2009) over the range 3.6 < 2 < 4.2 based on 
stacked SDSS quasar spectra. (Bottom) The blue dashed line is the 
model of Fig. 7 of Gnedin & Fan (2006) computed for the present 
cosmology. The black dot-dash line is the mfp computed from the 
LLS density of Choudhury (2009, Fig. 2(g)) with /3 = 1.3. 



only and there may be additional systematic errors. 
Rescaled to the present cosmology, the MHR value is 
25((1 -I- z)/A.5)~'^'^ Mpc which, while very similar in 
shape, is about a factor of two lower than the present 
results. The mean free path has a steeper slope than 
that of Faucher-Giguere and is reduced by a factor of 1.2 
at 2; = 5, increasing the emissivity required to ionize the 
IGM by the same factor. This increases the well-known 
discrepancy between the required and known sources of 



ionization at these high redshifts relative to the results 
of that paper. 

Fan et al. (2006) have made an estimate of the mfp 
based on the measured transmission and an assumed log 
normal density distribution function, following Miralda- 
Escude, Haehhnelt & Rees (2000). This procedure is 
very uncertain at the high redshifts where the transmis- 
sion is dominated by the extreme underdense tail of the 
distribution function, which is poorly described by the 
log normal function (Becker et al. 2007). Perhaps not 
surpisingly, given this and other major uncertainties in 
the interpretation, the Fan et al. (2006) estimates are 
steeper in shape and about a factor of 2.5 - 4 lower in 
normalization than the present measurements in the in- 
terval 4.5 < z < 6. 

The steep drop in the mfp with redshift to z = 6, 
if extrapolated to higher redshift, predicts a very small 
proper length value of around 5 Mpc at 2: = 6.5. If 
this is smaller than the size of the HII bubbles at the 
overlap epoch there could be important consequences for 
the evolution of the ionization rate (e.g. Furlanetto and 
Mcsingcr 2009). 

However, it is probable that as we move to earlier pe- 
riods the evolution of the geometry and ionization rates 
will combine to modify the evolution of the LLS, and 
simple extrapolations to higher redshifts might not be 
justified despite the smooth evolution from z = — 6. 
Full numerical simulations, capable of the resolution re- 
quired to describe the Lyman limit systems, are still ex- 
tremely challenging and subject to many uncertainties. 
In Figure 10 we compare the present observations with 
the simulation of Kohler & Gnedin (2007) and Gnedin & 
Fan (2006). These models are not a particularly good de- 
scription of the lower redshift mfp or n{z), presumably 
because of uncertainties in the evolution of the inter- 
galactic ionization parameter. However, they emphasize 
that the mfp could evolve smoothly to 2: = 6 and then 
change rapidly at higher redshifts. In this case the mod- 
els correspond to reionization just above z = 6 and the 
mean free path drops rapidly above this redshift. 

However, even the sign of the evolution in the mean 
free path is poorly determined at high redshift. Choud- 
hury & Ferrara (2005, 2006) and Choudhury (2009) have 
computed semi-analytic models based on the formalism 
of Miralda-Escudc (2003) with ionization histories cho- 
sen to match the high redshift constraints. These models 
provide a very good match to the 2 = 3 — 6 LLS density. 
However, as can be seen in Figure 10, they predict that 
the mfp flattens above 2 = 6 and the n{z) drops. In 
these models the bulk of the IGM volume is still ionized 
up to 2 = 8 and the change in the mfp corresponds to 
the transition away from popIII ionization at redshifts 
2 < 8. 

We conclude that it is probably wise not to extrapolate 
the mfp beyond the measured range. However, the shape 
and normalization of the mfp evolution at 2 < 6 may still 
provide useful tests of the higher redshift evolution since 
any model must satisfy the 2 = 6 boundary condition. 

This research was supported by the National Science 
Foundation under grant AST 06-07871. The authors 
acknowledge helpful conversations with Andrea Ferrara 
and Tirthankar Choudhury. 
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